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ABSTRACT 

We present high-precision timing observations spanning up to nine years for 37 millisecond pulsars moni¬ 
tored with the Green Bank and Arecibo radio telescopes as part of the North American Nanohertz Observatory 
for Gravitational Waves (NANOGrav) project. We describe the observational and instrumental setups used to 
collect the data, and methodology applied for calculating pulse times of arrival; these include novel methods 
for measuring instrumental offsets and characterizing low signal-to-noise ratio timing results. The time of ar¬ 
rival data are fit to a physical timing model for each source, including terms that characterize time-variable 
dispersion measure and frequency-dependent pulse shape evolution. In conjunction with the timing model fit, 
we have performed a Bayesian analysis of a parameterized timing noise model for each source, and detect evi¬ 
dence for excess low-frequency, or “red,” timing noise in 10 of the pulsars. For 5 of these cases this is likely due 
to interstellar medium propagation effects rather than intrisic spin variations. Subsequent papers in this series 
will present further analysis of this data set aimed at detecting or limiting the presence of nanohertz-frequency 
gravitational wave signals. 

Subject headings: Gravitational waves - Methods: data analysis - Pulsars: general 


1. INTRODUCTION 

The era of gravitational-wave astronomy is expected to be¬ 
gin within the next decade. It will be heralded by the first 
direct detection of gravitational waves as perturbations in the 
spacetime metric due to acceleration of massive objects. An¬ 
ticipated gravitational wave sources include merging systems 
of supermassive black hole binaries (SMBHBs) and neutron- 
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star binaries, as well as inflation-era relics (e.g., Grishchuk 
2005) and cosmic strings (e.g., Vilenkin & Shellard 1994). 
Several major experiments are underway in order to detect 
and characterize gravitational waves. One type of experiment 
is a pulsar timing array (PTA), in which a collection of radio 
pulsars is monitored, providing sensitivity to gravitational ra¬ 
diation in the nanohertz region of the spectrum (Hobbs et al. 
2010 ). 
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Millisecond pulsar rotation is very stable, and pulsars in 
relativistic binary systems have already been used to place the 
most stringent experimental constraints on strong-field grav¬ 
ity so far. The orbital decay observed in such binary sys¬ 
tems serves as compelling indirect evidence of gravitational 
radiation, as the observed orbital decay rates match the ex¬ 
pected rates due to loss of energy and angular momentum via 
emission of gravitational waves (Fonseca et al. 2014; Kramer 
et al. 2006; Weisberg et al. 2010). Sazhin (1978) and De- 
tweiler (1979) were the first to suggest that pulsar signals 
can be used to directly measure gravitational waves, particu¬ 
larly those produced by SMBHB mergers. Hellings & Downs 
(1983) extended this view and developed the notion that gravi¬ 
tational waves produce pulse time-of-arrival (TOA) shifts that 
are correlated amongst a set of pulsars. In principle, this al¬ 
lows the gravitational wave signal to be unambiguously sepa¬ 
rated from other astrophysical phenomena affecting measured 
TOAs - these would be specific to each pulsar, thus uncorre¬ 
lated between different objects. 

The North American Nanohertz Observatory for Gravita¬ 
tional Waves (NANOGrav) 1 is one of several PTA programs 
across the globe that collectively form the International Pul¬ 
sar Timing Array (IPTA; Hobbs et al. 2010). These collab¬ 
orations regularly monitor the most stable members of the 
millisecond pulsar (MSP) population distributed across the 
sky in order to achieve the highest sensitivity possible to¬ 
wards gravitational wave detection. While no direct detec¬ 
tion has been made so far, individual PTA programs have 
yielded upper limits on the amplitude A gw of the characteris¬ 
tic gravitational-wave strain h c due to a stochastic background 
in the nanohertz regime (van Haasteren et al. 2011; Demorest 
et al. 2013; Shannon et al. 2013; Lentati et al. 2015), where 

W)=A gw ( T A r ) (i) 

and / is the observed gravitational wave frequency. A num¬ 
ber of predictions for the expected strength of the SMBHB 
gravitational wave background have been made (e.g., Jaffe & 
Backer 2003; Sesana et al. 2008; Sesana 2013; McWilliams 
et al. 2014), with a = -2/3 and expected values for A gw rang¬ 
ing from ^5x 10 -16 to 4 x 10 _1 . Currently, the best pub¬ 
lished experimental limit is A gw < 2.7 x 10 -15 (Shannon et al. 
2013). In addition to measuring the gravitational wave back¬ 
ground, PTA measurements can be used to attempt to measure 
other gravitational wave signals, such as periodic gravitational 
waves from individual sources (Sesana et al. 2009; Arzouma¬ 
nian et al. 2014) and permanent deformations in spacetime, 
referred to as gravitational-wave “memory” (e.g., Madison 
et al. 2014; Wang et al. 2015). PTA data can also be used 
for ancillary studies of pulsars, binary systems, and the inter¬ 
stellar medium. 

In this study, we extend the data set analyzed by Demorest 
et al. (2013) and present high-precision timing observations 
of the updated NANOGrav PTA. The data comprise measure¬ 
ments of 37 MSPs and span nine years of observation. In 
Section 2, we provide information regarding the telescopes, 
methods and pulsar backends used for data collection. In Sec¬ 
tion 3, we describe the general procedure for data reduction, 
flux and polarization calibration, TOA determination and data 
excision. In Section 4, we outline the strategy used for gen¬ 
erating updated timing models for each pulsar and discuss the 

1 http://www.nanograv.org 


results of the model fitting. In Section 5, we describe the 
models used for characterizing noise in our timing data. In 
Section 6, we summarize the results and implications of this 
work. Raw and processed data products presented here are 
publicly available as of the date this work is published (§4). 

2. OBSERVATIONS 

This paper reports on observations of an array of millisec¬ 
ond pulsars made over a 9-year span from 2004 to 2013. Pul¬ 
sars were chosen for this project based on expectations of 
high time-of-arrival precision, reliable detection across a wide 
range of frequencies, and lack of unpredictable timing fluctu¬ 
ations from astrophysical effects (for example, no eclipsing 
binary pulsars have been included). The array initially in¬ 
cluded 15 pulsars, and it grew to 37 pulsars over the course of 
the project. The growth came from the discovery of new mil¬ 
lisecond pulsars and the advent of wide-band data acquisition 
systems, which allowed observation of some sources previ¬ 
ously deemed too weak or unreliable. The first five years of 
data on 17 of the pulsars were previously reported by Demor¬ 
est et al. (2013), however all data have been reprocessed as 
described in following sections. 

Pulsars with declinations in the range 0° < S < 39° were 
observed with the 305 m William E. Gordon Telescope of the 
Arecibo Observatory. Sources outside this declination range 
were observed with the 100 m Robert C. Byrd Green Bank 
Telescope (GBT) of the National Radio Astronomy Observa¬ 
tory. Two sources were observed at both telescopes, PSRs 
J1713+0747 andB1937+21. 

Table 1 summarizes the radio frequencies and data acquisi¬ 
tion systems used for this project; these are discussed in more 
detail below. Observation time spans for individual sources 
are listed in Table 2, and observation dates are displayed in 
Figure 1. 

Sources were observed at approximately monthly intervals 
through most of this program, with denser observations, ev¬ 
ery three weeks, in 2013. Scheduling of individual observing 
epochs varied depending on telescope operational consider¬ 
ations and sometimes deviated from regular cadences. Ob¬ 
servations at both telescopes were interrupted in 2007 due to 
telescope painting (Arecibo) and azimuth track refurbishment 
(GBT). 

Each pulsar was observed using radio receivers at two sep¬ 
arate frequencies throughout this program in order to mea¬ 
sure and remove frequency-dependent dispersive effects. Ex¬ 
ceptions were Arecibo observations of five sources that were 
made at a single frequency before 2009 or 2011 (depending 
on the source), and certain Arecibo observations during 2012, 
when technical issues impeded use of the 430 MHz receiver. 

At Arecibo, observations of a given pulsar using two re¬ 
ceivers were made in immediate succession within ^1 hour. 
At the GBT, observations using two receivers were typically 
separated by a few days due to the need for a physical re¬ 
ceiver change at that telescope. Observations without com¬ 
plementary data from the other receiver taken within 14 days 
were excluded from the data set. Exceptions to this rule were 
made for early Arecibo single-receiver observations of sev¬ 
eral sources mentioned above, and for wide-band 1400 MHz 
data, in which the wide frequency range of one receiver par¬ 
tially made up for the lack of data from a second receiver. The 
typical observation duration was about 25 minutes, with some 
variations over the course of the program. 

All receivers are sensitive to dual linear polarizations, with 
the exception of the Arecibo 430 MHz receiver, which mea- 
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Table 1 

Observing frequencies and bandwidths 


Telescope 

Receiver 


ASP/GASP 


PUPPI/GUPPI 


Data Span a 

Frequency 

Range b 

(MHz) 

Usable 

Bandwidth 0 

(MHz) 

Data Span a 

Frequency 

Range b 

(MHz) 

Usable 

Bandwidth 0 

(MHz) 

Arecibo 

327 

2005.0-2012.0 

315-339 

34 

2012.2-2013.8 

302-352 

50 

430 

2005.0-2012.3 

422-442 

20 

2012.2-2013.8 

421-445 

24 

L-wide 

2004.9-2012.3 

1380-1444 

64 

2012.2-2013.8 

1147-1765 

603 

S-wide 

2004.9-2012.6 

2316-2380 

64 

2012.2-2013.8 

1700-2404 d 

460 

GBT 

Rcvr_800 

2004.6-2011.0 

822-866 

64 

2010.2-2013.8 

722-919 

186 

Rcvrl_2 

2004.6-2010.8 

1386-1434 

48 

2010.2-2013.8 

1151-1885 

642 


a Dates of instrument use. Observation dates of individual pulsars vary; see Figure 1 . 

b Most common values; some observations differed. Some frequencies unusable due to radio frequency interfer¬ 
ence. 

c Nominal values after excluding narrow subbands with radio frequency interference. 
d Non-contiguous usable bands at 1700-1880 and 2050-2404 MHz. 


sures dual circular polarizations. Polarization cross-products 
were recorded so that full Stokes parameters could be re¬ 
covered. However, for the present work, we only use total- 
intensity measurements, obtained by summing the calibrated 
signals from pairs of orthogonal polarizations. 

Two sets of data acquisition systems were used. Early ob¬ 
servations (through 2012.3 at Arecibo and through 2011.0 at 
Green Bank) were recorded by the nearly-identical Astronom¬ 
ical Signal Processor (ASP) and Green Bank Astronomical 
Signal Processor (GASP) data acquisition systems (Demor- 
est 2007). Later observations (beginning 2012.2 at Arecibo 
and 2010.2 at Green Bank) were recorded using the nearly- 
identical Puerto Rican Ultimate Pulsar Processing Instrument 
(PUPPI) and the Green Bank Ultimate Pulsar Processing In¬ 
strument (GUPPI; DuPlain et al. 2008; Ford et al. 2010). Each 
of these systems digitized incoming baseband radio telescope 
voltage signals at the appropriate Nyquist rate, channelized 
them into subbands, performed real-time coherent dedisper¬ 
sion, calculated self- and cross-products to record full polar¬ 
ization information, and folded the data with the dynamically- 
calculated pulsar period using a pre-computed ephemeris. 
The end product of each instrument was folded pulse profiles 
(2048 bins) with self- and cross-products recorded over a se¬ 
ries of frequency channels and integrated over short subinter¬ 
vals over the course of an observation. 

The ASP and GASP systems processed up to 64 MHz of 
bandwidth, recording in contiguous 4 MHz subbands. Data 
were usually recorded in consecutive 60-second subintervals 
over the course of an observation. 

The PUPPI and GUPPI systems processed 100, 200, or 
800 MHz of bandwidth, depending on the mode of operation. 
In each case, data were recorded in contiguous subbands of 
width 1.5625 MHz. Data were usually recorded in consecu¬ 
tive 10-second subintervals, with 1-second subintervals used 
for some Arecibo 1400 MHz observations to aid interference 
excision. 

In some cases, particularly with PUPPI and GUPPI, band¬ 
width was limited by telescope receivers rather than the data 
acquisition instrument (Table 1). In post-processing, narrow- 
band radio frequency noise was removed and adjacent sub¬ 
bands were summed before arrival times were calculated (see 


§3.1). 

Each pulsar observation was preceded or followed by mea¬ 
surement of a pulsed noise signal using a setup identical to 
the pulsar observation in order to calibrate the signal levels. 
The pulsed noise signals themselves were calibrated in on- 
and off-source observations of unpolarized continuum radio 
sources on a monthly basis. 

For the timing analysis in this paper, only the polarization 
self-products were used. Data were summed in time and po¬ 
larization (§3.1). Simultaneous observations between ASP 
and PUPPI at Arecibo, and between GASP and GUPPI at 
the GBT, were used to measure the time offset between these 
pairs of instruments (Appendix A). 

3. CALIBRATION AND TOA DETERMINATION 

The results of the observations described in Section 2 are 
“raw” pulse profiles. This section describes the procedures 
employed to turn the raw profiles into usable pulse times of ar¬ 
rival (TOAs). These included: RFI excision, polarization and 
flux calibration, additional averaging in time and frequency, 
derivation of template profiles, and finally TOA determina¬ 
tion. All data processing operations described in this sec¬ 
tion were carried out using the psrchive software package 
(Hotan et al. 2004; van Straten et al. 2012). 2 These were 
organized into a processing pipeline via a set of scripts that 
are available online. 3 Overall, the calibration and processing 
strategy used here is similar to Demorest et al. (2013), and is 
based on standard methods for pulsar data analysis. 

3.1. Calibration and averaging 

The polarization calibration procedure used for these data 
sets was identical to that described by Demorest et al. (2013): 
A locally-generated broadband noise source is pulsed at 
25 Hz, split into two copies, and coupled in to the two polar¬ 
ization signal paths. Before each pulsar observation, a short 
(^1 minute) observation of the pulsed noise signal is recorded 
by the backend systems. This correlated noise source obser¬ 
vation is used to calibrate the two leading polarization terms: 

2 http://psrchive.sourceforge.net; PSRCHIVE version 
2015-01-15 b482 6eb was used for this work. 

3 http://github.com/demorest/nanopipe 
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J0023 + 0923 
J0030 + 0451 
J0340 + 4130 
J0613 — 0200 
J0645 + 5158 
J0931 —1902 
J1012 + 5307 
J1024 — 0719 
J1455 — 3330 
J1600 — 3053 
J1614 — 2230 
J1640 + 2224 
J1643 — 1224 

J1713 + 0747 

J1738 + 0333 
J1741 + 1351 
J1744 — 1134 
J1747 — 4036 
J1832 — 0836 
J1853 + 1303 
B1855 + 09 
J1903 + 0327 
J1909 — 3744 
J1910+ 1256 
J1918 — 0642 
J1923 + 2515 

B1937 + 21 

J1944 + 0907 
J1949 + 3106 
B1953+ 29 
J2010 —1323 
J2017 + 0603 
J2043 + 1711 
J2145 — 0750 
J2214 + 3000 
J2302 + 4442 
J2317 +1439 


2004 


2006 


2008 


2010 
Date [yr] 


2012 


2014 


Figure 1. Epochs of all observations in the data set. Marker type indicates data acquisition system: open circles are ASP or GASP; closed circles are PUPPI 
or GUPPI. Colors indicate radio frequency band, at either telescope: red is 327 MHz; orange is 430 MHz; green is 820 MHz; blue is 1.4 GHz; and purple is 
2.1 GHz. 
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differential gain and phase between the two hands of polar¬ 
ization. The noise source power is not assumed to be equal in 
each hand - its power in each polarization is measured sepa¬ 
rately at each observing epoch by observing the noise source 
while the telescope is pointed on and off a bright, unpolarized 
quasar (B1442+101 at Green Bank; J1413+1509 at Arecibo). 
For purposes of this paper, we used this calibration to balance 
the gains of orthogonal polarizations before summing to pro¬ 
duce total intensity profiles used for pulse timing. While we 
have not solved for a complete polarization calibration solu¬ 
tion here, the calibration data can also be used to create full- 
polarimetry profiles and to flux-calibrate the pulse profiles us¬ 
ing the known flux densities of the quasars used as calibration 
sources. 

Following calibration, excision of data corrupted by ra¬ 
dio frequency interference (RFI) was performed in two steps. 
First, a set of consistently bad frequency ranges for each tele¬ 
scope receiver was determined manually and was removed 
from all data sets. The interfering signals primarily respon¬ 
sible for these cuts are satellite transmissions near ~1.6 GHz 
and radar signals near ^1.2 GHz. This step resulted in re¬ 
moval of up to 15% of the full bandwidth recorded by GUPPI 
and PUPPI; the narrow-band ASP and GASP instruments 
were tuned to avoid these strong signals. The remaining us¬ 
able bandwidth for each receiver and data acquisition system 
is listed in Table 1 . Following this initial cut, remaining RFI 
was removed via the median filter algorithm in PSRCHIVE. 
In each 20-channel wide frequency window, the median off- 
pulse variation was computed, and any channels exceeding 4 
times this value were removed. Finally, prior to the final av¬ 
eraging described below, all profiles were normalized to have 
constant off-pulse variance. This step acted to down-weight 
any remaining corrupted data. 

As the TOA-determination procedure described in the next 
section begins to fail at very low signal-to-noise ratios (see 
Appendix B), it is advantageous to average as much data as 
possible into each profile before measuring a TOA. The fi¬ 
nal time and frequency resolution that should be used is ulti¬ 
mately limited by the need to resolve TOA shifts as a func¬ 
tion of time (for example, from orbital motion) or frequency 
(from profile frequency evolution or variable dispersion mea¬ 
sure). We chose to average profiles in time up to a maximum 
of 30 minutes or 2% of the pulsar’s binary period, whichever 
is shorter. Data from each observing session were divided 
equally - for example, a 40-minute observation would be 
averaged into two 20-minute sections. The shortest averag¬ 
ing was for PSR J0023+0923, which has an orbital period of 
200 minutes (see P/ ? column in Table 2). For frequency aver¬ 
aging, we adopted a slightly different strategy for each instru¬ 
ment. For ASP and GASP, no frequency averaging was done, 
and the final profile data remain at the instrumental resolution 
of 4 MHz. For GUPPI and PUPPI, the data were averaged 
to different final frequency resolutions depending on which 
receiver system was in use: 1.5625 MHz for 327 MHz and 
430 MHz data, 3.125 MHz for 0.82 GHz, and 12.5 MHz for 
all frequencies above 1 GHz. 

To summarize, the profile data set for any given observation 
consists of calibrated total-intensity profiles collected simul¬ 
taneously across many subbands (typically between 5 and 60) 
divided into one or more subintervals (typically 20-30 min¬ 
utes). 

3.2. Measuring times-of-arrival 


We calculated a pulse time-of-arrival (TOA) from each av¬ 
eraged profile resulting from the procedure described in the 
previous section. Thus any given observation results in a large 
number of TOAs, computed from data collected simultane¬ 
ously in different subbands. 

We calculated TOAs using the Fourier-domain algorithm 
of Taylor (1992) as implemented in the PSRCHIVE program 
pat. This method determined each TOA and its uncertainty 
via a least-squares fit for the pulse phase shift between an ob¬ 
served total-intensity pulse profile and an ideal template pro¬ 
file. Template determination was done using the same pro¬ 
cedure as Demorest et al. (2013). In short, for each pulsar 
and each receiver, we made a signal-to-noise-weighted sum 
of all GUPPI/PUPPI profile data. We then de-noised these 
profiles via wavelet decomposition and thresholding of the 
wavelet coefficients (as implemented in the PSRCHIVE pro¬ 
gram psrsmooth). The same template profiles were used 
to calculate TOAs from both GUPPI/PUPPI data and from 
GASP/ASP data. All templates were aligned so that phase 
zero occurs at the peak of the pulse profile. 

The pairs of data acquisition systems used at each tele¬ 
scope (GASP and GUPPI or ASP and PUPPI) had different 
signal path lengths and different internal latencies, which led 
to systematic TOA offsets. These must be measured and re¬ 
moved in order to avoid corrupting the pulsar timing results. 
This has typically been done in the past using the pulsar TOA 
measurements themselves - a time offset between two sys¬ 
tems is fit for, either as a term in the overall timing model fit 
(see Section 4), or separately using a subset of contempora¬ 
neous TOAs (e.g., Taylor & Weisberg 1989). More recently, 
Manchester et al. (2013) applied a method where a locally- 
generated timing signal was injected, measured and used to 
derive per-backend timing offsets. For the present work, we 
developed a new method that analyzes the noise in the pul¬ 
sar profiles collected simultaneously with a given pair of data 
acquisition systems. In simultaneous data the noise is corre¬ 
lated, and cross-correlating the pulse profile data from the two 
systems provided much higher-precision offset measurements 
than could be made from TOAs, where by design most of the 
noise was filtered out by the template-matching process. The 
results were offsets between GASP and GUPPI (also ASP and 
PUPPI), with typical value ~1 /zs and uncertainty ^5 ns, that 
were applied directly to the TOAs. Additional details are pre¬ 
sented in Appendix A. 

3.3. Editing time-of-arrival data sets 

After all TOAs were generated as described above, several 
cuts were made on the set of TOAs from each pulsar in order 
to arrive at the final set of data used in the analysis described 
in the next section. 

1. For observations where simultaneous data were 
recorded with both sets of instrumentation, the 
ASP/GASP TOAs were removed. 

2. In order to meaningfully determine a time-variable dis¬ 
persion measure, data from observing epochs with low 
fractional bandwidth, r'high/^iow <1.1, were removed. 
In practice, this criterion removed TOAs for any pul¬ 
sar that did not have, within any given 14-day window, 
either TOAs measured using two separate receivers or 
TOAs measured using one wide-band receiver with a 
wide-band data acquisition system. 
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3. As TOA measurement uncertainties become both un¬ 
derestimated and significantly non-Gaussian at low 
signal-to-noise ratio, TOAs from profiles with signal- 
to-noise ratio less than 8 were removed. See Ap¬ 
pendix B for further analysis of this cut. 

4. A small number of outlier TOAs were manually identi¬ 
fied and removed during the timing analysis. Typically 
these were due to low-S/N data that were just above 
the S/N cutoff or to RFI not excised by the algorithms 
described above. The TOAs removed by this process 
comprise ~1% of the full data set. 

All TOA data presented in this paper are publicly avail¬ 
able. 4 TOAs removed from the analysis as described above 
are included as supplementary files along with the main data 
set. The TOAs are given in TEMP02 format, with additional 
flags specifying relevant meta-data (e.g., backend, receiver, 
profile template file, etc). This data format can be read using 
both TEMPO 5 and TEMP02 6 pulsar timing analysis software. 
Clock correction data needed to reference the TOAs to the 
TT(BIPM) timescale (see Section 4.1) are distributed along 
with the TOA data, in the standard formats used by TEMPO 
and TEMP02. 

4. TIMING ANALYSIS 

We fit the TOA set for each pulsar to a timing model us¬ 
ing standard procedures as described by Lorimer & Kramer 
(2005), supplemented by novel methods to compensate for 
frequency-dependent pulse shape variations and to model 
noise in the TOA data sets. In this section, we start by sum¬ 
marizing our use of standard timing models, and then we de¬ 
scribe the new algorithm for handling frequency-dependent 
pulse shape variations. In the next section of the paper, we 
describe the noise model. 

4.1. Timing models 

All TOAs were initially measured using a local, topocen- 
tric time provided by hydrogen-maser clocks located at the 
observatories. Observatory clock corrections, determined by 
daily monitoring of the maser offsets compared to times de¬ 
termined using Global Positioning System (GPS) receivers, 
were used to transform the TOAs to Coordinated Universal 
Time (UTC). The times were further transformed to Terres¬ 
trial Time (TT) as published by the Bureau International des 
Poids et Mesures, TT(BIPM), after accounting for the ef¬ 
fects of the Earth’s varying rotation rate as published by the 
International Earth Rotation and Reference Systems Service 
(IERS). Finally, relativistic corrections are applied to convert 
the times to B ary centric Dynamical Time (TDB) 7 . Propaga¬ 
tion delays in the solar system, used to project the TOAs to 
the Solar-System Barycenter (SSB) are calculated using the 
JPL DE421 planetary ephemeris 8 , rotated into an ecliptic ref¬ 
erence frame using the 2010 IAU value of the obliquity of the 
ecliptic. 

4 http://data.nanograv.org 

5 http://tempo.sourceforge.net; TEMPO version 

2 014-11-2 0 7 6b837 5 was used for this work. 

6 http://tempo2.sourceforge.net; TEMP02 version 

2 014.11.1 was used for this work. 

7 See http: //www. iausofa.org/2 015_02 0 9_C/sofa/sofa_ 
ts_c.pdf for a detailed discussion of the various timescale transforma¬ 
tions. 

8 http://naif.jpl.nasa.gov/naif/ 


We derived timing models by enumerating all rotations of 
each pulsar and accounting for the various physical processes, 
discussed below, that can cause observed timing delays. This 
modelling was done in conjunction with a parameterized 
model for noise in the arrival-time data, described in detail 
in Section 5 and Appendix C. In effect, the noise model deter¬ 
mined a separate “weight” for each data subset, defined by the 
combination of receiver and backend system used, along with 
a measurement of temporally correlated “red” noise. We used 
the TEMPO and TEMP02 pulsar-timing analysis programs, 
making use of recently implemented generalized least-squares 
(GLS) fitting procedures which take into account correlations 
in the TOA noise when determining timing model parameter 
values and their uncertainties (e.g., Coles et al. 2011). 

The code-bases for tempo and tempo2 are different but 
not fully independent. With appropriate timing model options 
such that the two programs employed the same algorithms and 
conventions (using the same clock standards, employing the 
same solar system ephemeris, using the same obliquity of the 
ecliptic, excluding a solar wind model (see below), etc.) the 
fit results were nearly identical between these two programs. 
The vast majority of all fit parameters agreed to <10% of their 
1 -a uncertainties, and the ephemeris files we provide are able 
to be used in both programs. PSR J1713+0747 is an excep¬ 
tion, as its complicated timing model includes time-varying 
orbital geometry terms that are handled slightly differently by 
each program. In this paper we report only the TEMPO results 
for all the pulsars. 

The number of fit parameters in each timing model depends 
on the observed spin, astrometric, and environmental proper¬ 
ties of the given pulsar. We used ecliptic coordinates to fit 
for all astrometric parameters in order to reduce parameter 
covariances, and we fit for parallax for all pulsars, regard¬ 
less of whether the resulting fit value was physically mean¬ 
ingful (i.e., positive) or significant. We also fit for proper 
motion for all pulsars except for the two (PSRs J0931-1902 
and J1832-0836) which had observing timespans less than 
one year. The timing models contain fits for spin frequency 
and its first time-derivative, with higher-order spin noise, if 
present, being parametrized by the red noise model. We fit 
five Keplerian binary parameters for all binary pulsars using 
either the Damour & Deruelle (1985, “DD”) or Lange et al. 
(2001, “ELL1”) binary models. The former is a generally 
applicable, fully relativistic description of the pulsar’s orbit, 
while the latter is an alternate parameterization that improves 
numerical stability for very low-eccentricity orbits. We in¬ 
troduced post-Keplerian parameters (e.g., Damour & Taylor 
1992) when model accuracy was significantly improved as de¬ 
termined by an E-test significance value of 0.0027 (i.e., 3 -a 
significant). Information on the timing models, noise models 
and residual statistics is presented in Table 3. 

We incorporated timing model parameters to describe dis¬ 
persive delays in the TOAs from the time-variable integrated 
column of ionized gas between the observatory and each pul¬ 
sar. These delays are primarily due to the turbulent inter¬ 
stellar medium (ISM) but also include smaller contributions 
from the solar wind and the Earth’s ionosphere. The delays 
are characterized by the time-dependent dispersion measure 
(DM) for each pulsar which is directly proportional to the in¬ 
tegrated electron column density. The expected TOA delay 
for a broadband radio pulse is AfoM oc DMz/ -2 , where v is the 
observing frequency. 

We measured a value of DM at nearly every observing 
epoch for each pulsar, enumerated using the tempo/tempo2 
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Table 2 

Basic pulsar parameters and TOA statistics 


Source 

P 

dP/dt 

DM 

Pb 

Median scaled TOA uncertainty* 3 (ns) / Number of epochs 


Span 


(ms) 

(icr 20 ) 

(pc cm -3 ) 

(d) 

327 MHz 

430 MHz 

820 MHz 

1.4 GHz 

2.3 GHz 

(yr) 

J0023+0923 

3.05 

1.14 

14.3 

0.1 

- 

0.151 

22 

- 


0.179 

28 

- 


2.3 

J0030+0451 

4.87 

1.02 

4.3 

- 

- 

0.265 

53 

- 


0.261 

60 

- 


00 

bo 

J0340+4130 

3.30 

0.71 

49.6 

- 

- 

- 


0.782 

25 

1.595 

27 

- 


1.7 

J0613-0200 

3.06 

0.96 

38.8 

1.2 

- 

- 


0.116 

88 

0.450 

90 

- 


8.6 

J0645+5158 

8.85 

0.49 

18.2 

- 

- 

- 


0.269 

29 

0.985 

32 

- 


2.4 

J0931-1902 

4.64 

0.41 

41.5 

- 

- 

- 


0.804 

8 

1.554 

11 

- 


0.6 

J1012+5307 

5.26 

1.71 

9.0 

0.6 

- 

- 


0.377 

96 

0.518 

99 

- 


9.2 

J1024-0719 

5.16 

1.86 

6.5 

- 

- 

- 


0.582 

50 

0.846 

53 

- 


4.0 

J1455-3330 

7.99 

2.43 

13.6 

76.2 

- 

- 


0.868 

82 

1.766 

80 

- 


9.2 

J1600-3053 

3.60 

0.95 

52.3 

14.3 

- 

- 


0.267 

71 

0.202 

77 

- 


6.0 

J1614-2230 

3.15 

0.96 

34.5 

8.7 

- 

- 


0.336 

51 

0.424 

66 

- 


5.1 

J1640+2224 

3.16 

0.28 

18.5 

175.5 

- 

0.076 

61 

- 


0.082 

67 

- 


8.9 

J1643-1224 

4.62 

1.85 

62.4 

147.0 

- 

- 


0.301 

93 

0.481 

93 

- 


9.0 

J1713+0747 

4.57 

0.85 

16.0 

67.8 

- 

- 


0.100 

90 

0.050 

175 

0.025 

68 

00 

bo 

J1738+0333 

5.85 

2.41 

33.8 

0.4 

- 

- 


- 


0.316 

29 

0.301 

27 

4.0 

J1741+1351 

3.75 

3.02 

24.2 

16.3 

- 

0.155 

20 

- 


0.233 

42 

- 


4.2 

J1744-1134 

4.07 

0.89 

3.1 

- 

- 

- 


0.114 

89 

0.203 

91 

- 


9.2 

J1747-4036 

1.65 

1.32 

153.0 

- 

- 

- 


0.895 

22 

1.112 

25 

- 


1.7 

J1832-0836 

2.72 

0.87 

28.2 

- 

- 

- 


0.577 

11 

0.422 

10 

- 


0.6 

J1853+1303 

4.09 

0.87 

30.6 

115.7 

- 

0.369 

18 

- 


0.369 

50 

- 


5.6 

B1855+09 

5.36 

1.78 

13.3 

12.3 

- 

0.155 

74 

- 


0.148 

84 

- 


8.9 

J1903+0327 

2.15 

1.88 

297.6 

95.2 

- 

- 


- 


0.444 

35 

0.437 

35 

4.0 

J1909-3744 

2.95 

1.40 

10.4 

1.5 

- 

- 


0.041 

89 

0.102 

98 

- 


9.1 

J1910+1256 

4.98 

0.97 

38.1 

58.5 

- 

- 


- 


0.239 

75 

0.275 

37 

00 

bo 

J1918-0642 

7.65 

2.57 

26.6 

10.9 

- 

- 


0.317 

86 

0.547 

92 

- 


9.0 

J1923+2515 

3.79 

0.96 

18.9 

- 

- 

0.442 

17 

- 


0.535 

24 

- 


2.2 

B1937+21 

1.56 

10.51 

71.0 

- 

- 

- 


0.007 

93 

0.012 

154 

0.007 

47 

9.1 

J1944+0907 

5.19 

1.73 

24.3 

- 

- 

0.365 

24 

- 


0.403 

47 

- 


5.7 

J1949+3106 

13.14 

9.96 

164.1 

1.9 

- 

- 


- 


1.414 

22 

1.349 

15 

1.2 

B1953+29 

6.13 

2.97 

104.5 

117.3 

- 

0.475 

19 

- 


0.558 

51 

- 


7.2 

J2010-1323 

5.22 

0.48 

22.2 

- 

- 

- 


0.370 

53 

0.733 

55 

- 


4.1 

J2017+0603 

2.90 

0.80 

23.9 

2.2 

- 

0.237 

6 

- 


0.238 

27 

0.243 

20 

1.7 

J2043+1711 

2.38 

0.52 

20.7 

1.5 

- 

0.121 

24 

- 


0.170 

32 

- 


2.3 

J2145-0750 

16.05 

2.98 

9.0 

6.8 

- 

- 


0.213 

71 

0.535 

73 

- 


9.1 

J2214+3000 

3.12 

1.47 

22.6 

0.4 

- 

- 


- 


0.399 

26 

0.399 

22 

2.1 

J2302+4442 

5.19 

1.38 

13.7 

125.9 

- 

- 


1.018 

26 

1.592 

26 

- 


1.7 

J2317+1439 

3.45 

0.24 

21.9 

2.5 

0.070 75 

0.070 

74 

- 


0.202 

19 

- 


8.9 


Nominal scaling factor^ 7 (ASP/GASP) 

0.6 

0.4 


0.8 


0.8 


0.8 



Nominal scaling factor* (GUPPI/PUPPI) 

0.7 

0.5 


1.4 


2.5 


2.1 




a For this table, the original TOA uncertainties were scaled by their bandwidth-time product ( 10( f^ Hz 18 q 0 s ) to remove variation due to different instrument 

bandwidths and integration time. ' 

b TOA uncertainties can be rescaled to the nominal full instrumental bandwidth as listed in Table 1 by dividing by the scaling factors given here. 


parameter “DMX”. Since dual-receiver observations were 
sometimes separated by several days, we allowed a single 
constant DM value to apply to a window of up to 14 days 
of observations. These measured DM values include the ef¬ 
fects of ionospheric, interstellar, and solar wind dispersion 
(i.e., the solar wind DM model typically applied by default 
in TEMP0/TEMP02 was disabled for this analysis). The best- 
fit DMX values are shown in the timing summary figures for 
each pulsar below. The average DM value for each pulsar is 
highly covariant with profile shape evolution versus frequency 
(§4.2). However, DM variation can be easily distinguished 
from the constant-in-time profile shape terms. The DM error 
bars shown in the summary figures represent the uncertainty 
on the mean-subtracted values (DMX/- (DMX)), removing 
the large covariance affecting the mean DM. These uncertain¬ 


ties are determined via an appropriate linear transformation of 
the original post-fit parameter covariance matrix. 

Variations in DM over time are primarily attributed to an 
evolving view of the line-of-sight electron column due to the 
relative motion of the Earth, the pulsar, and the ISM (e.g. Ra- 
machandran et al. 2006). Many pulsars in our data set exhibit 
slow, long-term DM trends that are at least qualitatively in 
agreement with the expectation of turbulent electron density 
structure in the ISM. For several pulsars, annual variations in 
DM are also apparent. Those at low ecliptic latitude likely 
have a significant solar wind contribution to their DM - no¬ 
table examples are J1614-2230 and J2010-1323 with ecliptic 
latitudes of -1.2 and 6.5 degrees respectively. These pulsars 
show a sharp increase in DM at the times of year that they pass 
behind the Sun. Less sharp, more sinusoidal annual trends 
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are seen in some cases, notably J0613-0200 and J1643-1224. 
A possible explanation for this is line-of-sight motion due to 
Earth’s orbit projected onto an ISM density gradient, as dis¬ 
cussed by Keith et al. (2013). Using an independent data set, 
Keith et al. (2013) also detect significant annual DM modu¬ 
lation in several of the same sources where it is apparent in 
our results. Finally, isolated DM “events” such as the 2008- 
2009 dip shown by J1713+0747 are occasionally visible. This 
event is not yet well explained, but indicates the need for ad¬ 
ditional study. A more detailed astrophysical analysis of all 
DM results from this data set is currently underway. 

Regardless of the cause, since each pulsar’s DM is evolving 
with time, the fact that our dual-frequency measurements are 
separated by up to 14 days will result in a DM estimation 
error due to the DM being slightly different at the time of 
each of the two observations. This effect was studied in detail 
by Lam et al. (2015), who conclude that 14-day separations 
can result in timing errors on the order of ^50 ns, depending 
on the turbulent spectrum of electron density fluctuations in 
the ISM. In practice this effect will be significantly smaller as 
the majority of our paired observations are separated by <3 
days (e.g., Lam et al. 2015, Figure 4). 

4.2. Compensation for frequency evolution of pulsar profiles 

For each pulsar, we calculated TOAs for all frequency chan¬ 
nels recorded with a given receiver using a single standard 
template profile. Because pulse shapes vary with frequency, 
this produced small systematic frequency-dependent pertur¬ 
bations in the TOAs in addition to the u~ 2 offsets due to dis¬ 
persion. An example of such common behavior in our data set 
is in Figure 2. In previous work, we compensated for these 
frequency-dependent perturbations for each pulsar by fitting 
arbitrary time offsets to each spectral channel, thus adding a 
large number of free parameters to the timing model (Demor- 
est et al. 2013). 

For this analysis, we developed a heuristic approach to re¬ 
move as much of the frequency-dependent (“FD”) bias from 
our timing residuals as possible by incorporating an additional 
timing delay Af FD to all timing models, where 

a ' ro= i>i° 8 (iSfe) <2> 

i= 1 v 7 

and the coefficients c t are fit parameters in the timing models. 
For any given pulsar, the number of terms needed was deter¬ 
mined by an F-test significance value of 0.0027, the same cri¬ 
terion used for other timing model parameters. The number of 
parameters used ranged from 0 to 5 (Table 3), and reflects the 
degree to which profile evolution is important for each pulsar. 
An illustrative example of the application of our FD model for 
bias removal is shown in Figure 2 (see also discussion of this 
approach by Zhu et al. 2015). 

This algorithm significantly reduced the number of free pa¬ 
rameters in the timing model compared to our previous work, 
but it remains an ad hoc procedure - it is applied only after 
TOAs are calculated and does not directly utilize the addi¬ 
tional pulse shape information available in the profile data. In 
future work we plan to explore additional approaches to solv¬ 
ing this problem, including TOAs derived from a broadband 
profile-fitting approach that directly incorporates profile evo¬ 
lution versus frequency (Pennucci et al. 2014). 

4.3. Timing summary 



600 1000 1400 1800 2200 


Frequency (MHz) 

Figure 2. Average timing residual versus radio frequency for 
PSR J1713+0747. Upper panel: Non-dispersive frequency dependent 
(FD) residuals with all other timing model parameters held fixed at their 
best-fit values. The dashed line indicates the best-fit FD model from Eqn. 2, 
offset for plot clarity. Middle panel: Residuals when FD model is not 
included in the fit. Lower panel: final residuals when FD model is included 
in the fit. Note the different y-axis scales in each panel. See also a similar 
presentation of these data by Zhu et al. (2015). 

A basic summary of the timing model fit results, includ¬ 
ing number of TOAs, number of fit parameters, noise model 
results, and basic statistics of the residuals is presented in Ta¬ 
ble 3. The full set of best-fit timing model parameter values, 
and their associated uncertainties, are publicly available and 
are distributed along with the TOA data presented here. The 
models are provided in the standard “par file” format under¬ 
stood by TEMPO and TEMP02. These results potentially con¬ 
tain a signficant amount of astrophysical information about 
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these pulsars, their orbits and binary companions, and the ISM 
properties along their line-of-sight, however we have post¬ 
poned a detailed interpretation of this to future work (see Sec¬ 
tion 6 for a partial list). It should be noted that while the 
models presented here provide an accurate description of the 
time of arrival data for the purposes of gravitational wave de¬ 
tection, certain model parameters (in particular, parallax and 
Shapiro delay) may require a more sophisticated uncertainty 
analysis before astrophysically meaningful conclusions can 
be drawn from their values. 

5. NOISE CHARACTERIZATION 
5.1. Noise Model 

The noise model used in the analysis is a parameterized one 
that includes the effects of several noise sources that produce 
different correlations of TOAs obtained in non-overlapping 
time blocks and frequency channels. For instance, the tem¬ 
plate matching errors due to radiometer noise are uncorrelated 
in both time and frequency, but pulse-jitter noise (Cordes & 
Shannon 2010) appears to affect all TOAs obtained simulta¬ 
neously in different frequency channels. Correlated timing 
noise with a red power spectrum occurs to varying degree 
in different pulsars. Spin noise is achromatic and is much 
smaller in MSPs compared to objects with stronger magnetic 
fields and longer spin periods. Chromatic red noise due to 
propagation through intervening plasmas (ISM, interplanetary 
medium, and ionosphere) may also be present if dispersive 
delays are not removed perfectly or if scattering and refrac¬ 
tion effects contribute significantly. Jitter noise appears to be 
highly correlated across hundreds of MHz for those MSPs that 
have been analyzed in detail (Shannon et al. 2014; Dolch et al. 
2014). An approach to noise modeling has been discussed 
extensively in van Haasteren & Levin (2013); Ellis (2013); 
van Haasteren & Vallisneri (2014, 2015); Arzoumanian et al. 
(2014); Ellis (2014). For more details about the specific noise 
model and implementation for the current study see Appendix 
C. Here we will summarize the parameterization used for this 
data release. 

Our model for noise starts with the measurement error on 
each TOA, a^ k , determined by the template-matching TOA 
calculation algorithm; here i is the TOA number and the sub¬ 
script k denotes the backend/receiver system. Because such 
measurement errors may be underestimated, we allow for 
them to be increased by systematic quadrature and scaling 
factors, determined on a system-by-system basis, 

<7 uk ^E k {c7l k + Ql) l/ \ (3) 

where E k and Q k are the EFAC and EQUAD parameters used 
in the TEMPO/TEMP02 timing code. 

In addition to the template-fitting errors, we also include 
TOA errors that are uncorrelated in time but completely cor¬ 
related between TOAs obtained at different frequencies mea¬ 
sured simultaneously, termed “short-term correlated noise” in 
Appendix C. The strength of this process is characterized by 
the ECORR parameter. This term could include true pulse 
phase jitter (Cordes & Shannon 2010), known to be present in 
some pulsars, but can also include other similarly-correlated 
components. Lastly, we model the red noise as a stationary 
Gaussian process that is parameterized by a power spectrum 
of the form / f \ 7red 

P(f)=AL(jr) , ( 4 ) 

where A re d is the amplitude of the red noise in /is yr 1 / 2 , 7 re d is 


the spectral index, and / yr = lyr -1 . 

These noise parameters are included in a joint likelihood 
that contains all timing model parameters. For the purposes 
of this paper, we analytically marginalize over the linear tim¬ 
ing model parameters and explore the space of noise pa¬ 
rameters via Markov Chain Monte-Carlo (MCMC). We then 
use the MCMC results to determine the maximum likelihood 
noise parameters, which are subsequently used as inputs to 
TEMPO/TEMP02 GLS fitting routines. For each pulsar we al¬ 
ways include the EFAC, EQUAD, and ECORR parameters 
for each backend/receiver system. We only include red noise 
when it is preferred by the data. The red noise model selec¬ 
tion is performed with MultiNest (Feroz et al. 2009) using a 
Bayes factor threshold of 100 to determine whether red noise 
is included in the final model. The applicable red noise ampli¬ 
tudes, spectral indices, and Bayes factors are shown in Table 
3. 

5.2. Noise Analysis 

Here we discuss some of the major features that our noise 
model reveals. As mentioned above, we only include red 
noise in the noise model when the data favor its inclusion. 
In our analysis, ten pulsars meet this criterion and will now 
be discussed further. Intrinsic pulsar spin noise and its ef¬ 
fects have been studied in (Blandford et al. 1984; Cordes 
& Downs 1985; Arzoumanian et al. 1994) and (Shannon & 
Cordes 2010, hereafter SC10). Using a sample of both canon¬ 
ical pulsars (CPs) and MSPs, SC 10 parameterize the post-fit 
(after quadratic subtraction) timing noise rms as 

0TN,2 = C2^ Q! |z/-l5| /3 T y 7 , (5) 

where v, z>_i 5 , and T yr are the spin frequency, spin frequency 
derivative in units of 10 -15 s -2 , and time span of the data set 
in years. The best fit values of the free parameters were found 
to be ln(C 2 ) = 1.6±0.4, a =-1.4±0.1,/3 = 1.1 ± 0 . 1 , and 7 = 
2.0 ± 0.2. A fifth parameter, S , was used to take into account 
the empirical scatter about the mean relation in Eq. 5 and was 
estimated to be S = 1.6 ± 0.1 in In cttn, 2 - First we note that the 
best fit value of 7 in Eq. (5) corresponds to a power spectral 
density index of 7 re d = -( 27 +1) = -5 ± 0.4 in Eq. (4). We can 
estimate <ttn ,2 from our noise model by 



( 6 ) 


where the lower integration limit of 1 /T serves as a filter for 
quadratic subtraction. Furthermore, we can produce a distri¬ 
bution of <ttn ,2 by evaluating Eq. ( 6 ) for all values of A re 4 and 
7 re d from our MCMC analysis. This will represent our un¬ 
certainty in the red noise variance by incorporating the full 
posterior distributions of red noise parameters as opposed to 
just the maximum likelihood values. 

In essence, SC 10 make two predictions for intrinsic pul¬ 
sar timing spin noise: (/) the red noise spectral index is steep 
with a value ^ -5, and (ii) the red noise rms follows Eq. (5) 
to within a factor of exp (±5). In Figure 3 we show the max¬ 
imum a-posteriori value and 68 % credible interval for the red 
noise spectral index, 7 red , for all pulsars that display signifi¬ 
cant red noise. 
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Table 3 

Summary of timing model fits 


Source 

Number 

ofTOAs 

Number of Fit Parameters' 3 

S A B DM FD J 

RMS fc (/xs) 
Full White 

Red Noise' 

Aed 7red 

c 

logio B 

Figure 

Number 

J0023+0923 

4598 

3 

5 

5 

28 

1 

1 

0.320 

- 

- 

- 

1.40 

7 

J0030+0451 

2468 

3 

5 

0 

60 

0 

1 

0.723 

0.212 

0.014 

—4.8 

4.99 

8 

J0340+4130 

3008 

3 

5 

0 

27 

1 

1 

0.385 

- 

- 

- 

0.01 

9 

J0613-0200 

7651 

3 

5 

7 

90 

2 

1 

0.592 

0.165 

0.093 

-2.9 

3.71 

10 

J0645+5158 

2896 

3 

5 

0 

33 

2 

1 

0.052 

- 

- 

- 

-0.08 

11 

J0931-1902 

719 

3 

3 

0 

11 

0 

1 

0.381 

- 

- 

- 

-0.07 

12 

J1012+5307 

11995 

3 

5 

5 

95 

1 

1 

1.197 

0.355 

0.669 

-1.0 

14.31 

13 

J1024-0719 

5073 

3 

5 

0 

53 

2 

1 

0.280 

- 

- 

- 

0.08 

14 

J1455-3330 

5122 

3 

5 

6 

81 

1 

1 

0.694 

- 

- 

- 

0.05 

15 

J1600-3053 

8174 

3 

5 

8 

74 

2 

1 

0.197 

- 

- 

- 

-0.01 

16 

J1614-2230 

7517 

3 

5 

7 

54 

1 

1 

0.189 

- 

- 

- 

0.02 

17 

J1640+2224 

2565 

3 

5 

9 

65 

2 

1 

0.158 

- 

- 

- 

-0.03 

18 

J1643-1224 

7119 

3 

5 

6 

91 

2 

1 

2.057 

0.331 

1.231* 

-1.7 

18.33 

19 

J1713+0747 

15830 

3 

5 

8 

106 

4 

3 

0.116 

- 

- 

- 

0.01 

20 

J1738+0333 

2711 

3 

5 

5 

28 

1 

1 

0.308 

- 

- 

- 

0.01 

21 

J1741+1351 

1600 

3 

5 

8 

27 

0 

1 

0.103 

- 

- 

- 

-0.02 

22 

J1744-1134 

9020 

3 

5 

0 

88 

2 

1 

0.334 

- 

- 

- 

0.25 

23 

J1747-4036 

2778 

3 

5 

0 

25 

1 

1 

0.531 

- 

- 

- 

0.12 

24 

J1832-0836 

1136 

3 

3 

0 

10 

0 

1 

0.121 

- 

- 

- 

-0.04 

25 

J1853+1303 

1411 

3 

5 

6 

26 

0 

1 

0.235 

- 

- 

- 

-0.02 

26 

B1855+09 

4071 

3 

5 

7 

72 

3 

1 

1.339 

0.505 

0.017 

-4.9 

2.87 

27 

J1903+0327 

1887 

3 

5 

8 

36 

2 

1 

1.949 

0.327 

0.851* 

-2.5 

2.87 

28 

J1909-3744 

10697 

3 

5 

8 

88 

1 

1 

0.079 

- 

- 

- 

0.72 

29 

J1910+1256 

2690 

3 

5 

6 

45 

1 

1 

1.449 

0.587 

0.801“* 

-1.9 

5.39 

30 

J1918-0642 

10035 

3 

5 

7 

87 

3 

1 

0.340 

- 

- 

- 

-0.02 

31 

J1923+2515 

939 

3 

5 

0 

24 

1 

1 

0.266 

- 

- 

- 

-0.06 

32 

B1937+21 

9966 

3 

5 

0 

102 

5 

3 

1.549 

0.104 

0.197 

-2.4 

96.48 

33 

J1944+0907 

1724 

3 

5 

0 

28 

2 

1 

2.442 

0.331 

0.860* 

-2.8 

2.35 

34 

J1949+3106 

1416 

3 

5 

7 

16 

0 

1 

0.647 

- 

- 

- 

0.03 

35 

B1953+29 

1329 

3 

5 

6 

24 

2 

1 

4.149 

0.531 

0.015“* 

-6.7 

2.14 

36 

J2010-1323 

8068 

3 

5 

0 

55 

1 

1 

0.312 

- 

- 

- 

0.08 

37 

J2017+0603 

1589 

3 

5 

7 

24 

0 

2 

0.073 

- 

- 

- 

0.01 

38 

J2043+1711 

1394 

3 

5 

7 

23 

1 

1 

0.108 

- 

- 

- 

0.02 

39 

J2145-0750 

7369 

3 

5 

6 

73 

2 

1 

0.371 

- 

- 

- 

0.20 

40 

J2214+3000 

2624 

3 

5 

5 

25 

1 

1 

0.319 

- 

- 

- 

0.06 

41 

J2302+4442 

3044 

3 

5 

6 

27 

1 

1 

0.708 

- 

- 

- 

0.47 

42 

J2317+1439 

2650 

3 

5 

8 

68 

3 

2 

0.267 

- 

- 

- 

0.04 

43 


a Fit parameters: S=spin; B=binary; A=astrometry; DM=dispersion measure; FD=frequency dependence; J=jump 

b Weighted root-mean-square of epoch-averaged post-fit timing residuals. For sources with red noise, the “Full” RMS value includes the red noise contribution, 

while the “White” RMS does not. 

c Red noise parameters: A rec j = amplitude of red noise spectrum at /= 1 yr -1 measured in /is yr 1 / 2 ; 7 rec i = spectral index; B = Bayes factor. See Eqn. 4 and 

Appendix C for details. 

d For these sources, the detected red noise is likely due to unmodeled interstellar medium propagation effects rather than intrinsic spin noise; see text for details. 


We see that our noise analysis yields a much more shallow 
spectral index in general than the predicted value of SC 10. In 
fact, of the 10 pulsars that display red noise, only 3 (PSRs 
J0030+0451, B1953+29, and B1855+09) have spectral in¬ 
dices consistent with -5, the others are consistent with ^ -2. 
If we assume that this red noise is due to a random walk in 
one of the quadratic spin down parameters, then our analy¬ 
sis suggests a random walk in the pulsar phase 9 . However, 
it is more likely that in many cases (pulsars marked with an 
asterisk in Figure 3) this behavior is due to un-modeled ISM 
effects as we will discuss. In the right panel of Figure 3 we 
see that our measurements of <ttn ,2 are close to one-sigma 

9 As stated in SC 10, random walks in the pulsar phase, period and pe¬ 
riod derivative lead to underlying power spectral indices of -2, -4, and -6, 
respectively. 


consistent with the predictions of SC 10 with the exception of 
PSR B1937+21 which exhibits much weaker red noise than 
predicted. The gray points show that the 95% upper limits on 
<jtn ,2 are not consistent for some pulsars. Overall we can state 
with confidence that our noise analysis is inconsistent with the 
predictions of SC 10 both for the spectral index and the overall 
red noise rms. To explore this more closely, we will now look 
into each pulsar in more detail. 

PSRs J0030+0451 and B1855+09 are consistent with the 
spin noise predictions of a steep red noise process. From in¬ 
spection of Figures 8 and 27 we see that both pulsars are timed 
for the full nine years and have dual frequency data and DM (t) 
corrections for all observing epochs. Furthermore, each set of 
residuals displays a cubic low frequency term that is charac¬ 
teristic of the predicted steep red process. This appears to be 
the first evidence of red noise in these pulsars as they have 
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Figure 3. Left panel: Maximum a-posteriori value and 68% credible interval on the red noise spectral index, 7 rec j for all pulsars that display significant red noise. 
The dashed and dotted black lines represent the mean and one-sigma predictions on the spectral index from SC 10. The points with square markers are the spectral 
index values presented in Lentati et al. (2015). The pulsars marked with an asterisk indicate those pulsars for which the red noise is likely due to unmodeled 
interstellar medium propagation effects, rather than intrinsic spin noise. Right panel: Measured value for red noise RMS for all red pulsars (see text for details of 
this calculation) vs. predictions for SC10. The uncertainties are the one-sigma and 68% credible region for the predicted and measured values, respectively. The 
gray points are the 95% upper limits from the predicted and measured values of cjtn ,2 for all other pulsars. The color scheme is the same as the left panel. We do 
not include the pulsars marked with an asterisk in this figure. 


white residuals for both five and six year datasets presented 
from NANOGrav (Demorest et al. 2013) and the Parkes Pul¬ 
sar Timing Array (PPTA Manchester et al. 2013). 

PSR B1953+29 is also consistent with the spin noise pre¬ 
dictions of a steep red noise process. However, as shown in 
Figure 36, this pulsar lacks dual-frequency data early in the 
timing campaign, which is likely a strong contributor to the 
measure red noise. 

Both PSRs J1643-1224 and J1910+1256 were identified as 
displaying strong evidence of red noise in Demorest et al. 
(2013), and PSR J1643-1224 had a significant v in Manch¬ 
ester et al. (2013). In the case of PSR J1643-1224 the shal¬ 
low red noise process may be due to uncorrected ISM effects 
that include scattering and refraction (e.g., Rickett 1990; Fos¬ 
ter & Cordes 1990; Cordes & Shannon 2010) - as can be seen 
in Figure 19, there is a clear dependence of the noise on ra¬ 
dio frequency. The red noise present in the residuals of PSR 
J1910+1256 is likely caused by DM variations due to the fact 
that we only have single-frequency observations for the first 
four years of the data set. DM variations for a Kolmogorov 
spectrum would give an / -8 / 3 spectrum of TOA variations but 
this can be altered by linear changes in DM from large-scale 
structures or changes in pulsar distances from their line-of- 
sight motions. This is further indicated by inspection of Fig¬ 
ure 30 where the timing residuals appear relatively white after 
dual frequency observations had begun. 

PSR B1937+21 displays the strongest red noise in our sam¬ 
ple, consistent with previous work that shows a large amount 
of red noise (e.g., Shannon & Cordes 2010; Shannon et al. 
2013). Unlike previous work, which indicates a steep, red 
spectrum our analysis shows a shallower spectrum. Although 
it has a large cubic trend, the shallow red noise spectral in¬ 
dex measurement indicates that there are still high frequency 
trends in the data. In fact, in Figure 33 we can see that the 
red noise seems to track the DM changes around 2011.5 - 
2012.5. This feature, along with the large DM (71 pc cm -3 ), 
suggest that unmodeled ISM effects may contribute to the ob¬ 
served red noise, particularly at high frequencies, resulting in 
a lower measured spectral index. An additional explanation 


for the much shallower spectral index is that we are analyzing 
only the most recent nine years of data on this pulsar whereas 
previous analyses have used a much longer time span (SC 10 
uses up to 24 years of data) not encompassing this new data. 
This indicates that the noise could be non-stationary in na¬ 
ture. PSR J1903+0327 has a similar feature around the same 
time in which a large drop in the DM coincides with a peak 
in the red noise. Once again, this effect, in combination with 
the very large DM (297.6 pc cm -3 ) indicate ISM effects as 
opposed to intrinsic instability as the cause of the red noise. 

The measured red noise in PSR J1944+0907 is, likely due 
to unmodeled DM or scattering/refraction effects since there 
is only single-frequency data for a large portion of the data 
set. 

As shown in Figure 3, PSRs J0613-0200 and J1012+5307 
also display low spectral-index red noise; however, although 
there are clear high frequency fluctuations in the residuals 
(Figs 10, 13) there are no obvious radio frequency dependent 
features present. Therefore, it is difficult to assess the cause of 
this measured red noise for these pulsars. Nonetheless, these 
results are very consistent with Lentati et al. (2015) (square 
marker in Figure 3) where a nearly identical noise model was 
used for these two pulsars observed with EPTA telescopes. 
This at the very least can rule out any instrumental effects. 

It was also pointed out in Manchester et al. (2013) that 
PSR J1909-3744 displayed some evidence of red noise for 
the PPTA and EPTA data sets. While our analysis of PSR 
J1909-3744 did note find sufficient red noise to classify it as 
a detection (Bayes factor greater than 100; Table 3), the pos¬ 
terior probability distributions for this pulsar hint at the pres¬ 
ence of weak red noise, again with a shallow spectral index. 
This is interesting in that PSR J1909-3744 has very good tim¬ 
ing precision and is ideal for GW detection prospects. Future 
longer data sets will test whether this pulsar truly displays red 
noise. 

Finally, we point out PSRs J1600-3053 and J1747-4036. 
While not displaying evidence of red noise, we do see non¬ 
white features in the residuals that are radio frequency depen¬ 
dent. Time varying DM corrections were included for the full 
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range of both of these datasets indicating that while the noise 
is radio frequency dependent, it is not likely a DM effect. 

The noise model used for this data release provides a good 
fit to the data for most of the pulsars in our data set. How¬ 
ever, the model does not accommodate time-variable, chro¬ 
matic phenomena other than DM variations with their u~ 2 
dependence on radio frequency. Such phenomena may in¬ 
clude frequency dependent dispersion measures (Cordes et al. 
2015) or scattering, and presently such phenomena are im¬ 
precisely absorbed by the red-noise and short-term correlated 
noise models. As the timespan of wide-band millisecond pul¬ 
sar data sets grows, it will become practical to incorporate 
such phenomena into the noise model. 

6 . SUMMARY AND CONCLUSIONS 

In summary, we have obtained, reduced, analyzed, and 
made public pulse times of arrival for 37 millisecond pulsars, 
using two telescopes, over a time span of up to nine years. 
A major upgrade in backend instrumentation occurred mid¬ 
way through the data set; we developed a novel method for 
measuring the instrumental offset between these systems that 
removes the need to fit this effect using the TOA data. We 
have continued to develop and refine methods for characteriz¬ 
ing time-variable dispersion measure and frequency-depenent 
pulse shape evolution, while fitting phyiscal timing models 
to these data. A significant new development is the parame¬ 
terized noise model presented in this paper, and its inclusion 
in the timing model fit via a generalized least squares proce¬ 
dure. Our noise modelling has indicated the presence of time- 
correlated, or red, noise in 10 of these pulsars; we suspect a 
combination of propagation effects in the interstellar medium 
and intrinsic spin noise both contribute to these detections, 
with levels of each varying on a case-by-case basis. 

The primary scientific motivation for this project is to detect 
or limit the presence of nanohertz-frequency gravitational ra¬ 
diation by looking for correlated timing fluctuations amongst 
this set of pulsars. While the analysis presented here deals 
with each pulsar separately, subsequent papers in this series 
will perform correlation analyses to look for the effect of 
different gravitational wave signal types. These include the 
stochastic background from supermassive black hole bina¬ 
ries and/or cosmic strings; continuous-wave emission from 
individual binary systems; and gravitational wave bursts with 
memory following merger events. 

In addition to the gravitational wave analyses just men¬ 
tioned, a number of additional topics are planned to be ad¬ 
dressed in future papers, including: Detailed investigation of 
pulse jitter and other sources of noise in these data; measure¬ 
ment of orbital parameters, pulsar and companion masses, 
and relativistic orbital effects; the effect of scattering on the 
timing results; pulsar astrometry, distance measurements and 
kinematics; analysis of the polarization properties of the pulse 
profiles; flux densities and population analysis; and reanaly¬ 
sis of these data using wide-band timing methodologies. Fu¬ 
ture improvements to the data set presented here include on¬ 
going increase in the number of pulsars measured, increased 
cadence on several of the best pulsars in the set, and inclusion 
of archival Arecibo data covering up to 20 years total times¬ 
pan. 

Author contributions. An alphabetical-order author list was 
used for this paper in recognition of the fact that a large, 
decade-timescale project such as this is necessarily the result 


of the work of many people. All authors contributed to the 
activities of the NANOGrav collaboration leading to the work 
presented here, and reviewed the manuscript text and figures 
prior to the paper’s submission. Additional specific contribu¬ 
tions to this paper are as follows. ZA, KC, PBD, TD, RDF, EF, 
MEG, GJ, MJ, MTL, LL, MAM, DJN, TTP, SMR, IHS, KS, 
JKS, and WWZ made observations for this project and devel¬ 
oped timing models (§4). PBD co-developed the GUPPI and 
PUPPI instruments, wrote observing proposals, coordinated 
GBT observations, performed the calibration and TOA gen¬ 
eration, developed and implemented the offset-measurement 
analysis, developed and implemented the FD profile evolution 
model, analyzed TOA accuracy in the low-signal-to-noise- 
regime, compiled the data files for public release, coordinated 
the writing of the paper, and contributed substantially to the 
text, figures, and tables. JAE developed and implemented 
the noise model, used it to analyze the present data set, an¬ 
alyzed the red-noise results, co-developed the daily-average 
residual statistic and implemented it for the residual figures, 
wrote substantial amounts of the text, and contributed tables 
and figures. DJN co-authored observing proposals, coordi¬ 
nated the NANOGrav Timing Working Group, wrote portions 
of the text, and contributed a figure and table. TTP undertook 
careful analysis of the residuals leading to the implementation 
of the signal-to-noise cutoff. IHS coordinated the Arecibo ob¬ 
servations and wrote portions of the text. SMR co-developed 
the GUPPI and PUPPI instruments and wrote portions of the 
text. RvH contributed to the noise model and analysis and 
verified it using independent code and co-developed the daily- 
average residual statistic. EF wrote portions of the text. ZA, 
JMC, TD, and others edited the manuscript text and/or pro¬ 
vided significant input about its content. 

The NANOGrav project receives support from National 
Science Foundation (NSF) PIRE program award number 
0968296 and NSF Physics Frontier Center award number 
1430284. NANOGrav research at UBC is supported by 
an NSERC Discovery Grant and Discovery Accelerator 
Supplement and the Canadian Institute for Advanced Re¬ 
search. Part of this research was carried out at the Jet 
Propulsion Laboratory, California Institute of Technology, 
under a contract with the National Aeronautics and Space 
Administration. TTP is a student at the National Radio 
Astronomy Observatory. The National Radio Astronomy 
Observatory is a facility of the NSF operated under co¬ 
operative agreement by Associated Universities, Inc. The 
Arecibo Observatory is operated by SRI International under 
a cooperative agreement with the NSF (AST-1100968), and 
in alliance with Ana G. Mendez-Universidad Metropolitana, 
and the Universities Space Research Association. Some 
computational work was performed on the Nemo cluster 
at UWM supported by NSF grant No. 0923409. JAE 
and RvH admowledge support by NASA through Einstein 
Fellowship grants PF4-150120 and PF3-140116, respectively. 

APPENDIX 

A. TIMING OFFSET DETERMINATION 

A common problem in long-term pulsar timing studies is 
to connect the timing results between multiple generations of 
instrumentation at a given telescope. Single-dish pulsar tim¬ 
ing data is generally timestamped either at the point where 
the radio signal is digitized, or somewhat later at the output 
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of a filterbank, when the data are received by software sys¬ 
tems. Prior to this point, the signal accumulates additional 
latency as it passes through various telescope electronics sys¬ 
tems. These can include analog cable delays (tens to hundreds 
of ns), transmission over long fiber optic links (tens of /is), 
and filter latencies (up to a few /is). Since only variations in 
pulse phase - not its absolute value - are physically meaning¬ 
ful in timing analyses, the presence of time offsets like these 
is not a problem, as long as they are constant. However, when 
a new backend instrument is added, the delay values (e.g., ca¬ 
bling, filters) typically differ from the previous version. As 
described in §3, this time offset must be accounted for in or¬ 
der to measure long-timescale effects. 

A typical approach used in most past timing analyses is 
to measure offsets between instruments using pulsar TOAs 
(e.g., Taylor & Weisberg 1989). An arbitrary offset can 
be included as an additional term (known as a “JUMP” in 
TEMP0/TEMP02) in a timing model fit (see §4). Although 
this is a straightforward approach, it has several drawbacks: 
First, in the presence of unmodelled red noise, this can in¬ 
troduce systematic biases in other parameters (see for exam¬ 
ple discussion in Coles et al. 2011); this can be mitigated via 
improved noise modelling as described in §5 and references 
therein. Next, even if the noise is properly modelled, the offset 
will be covariant with other long-term effects, most critically 
with low-frequency gravitational waves, reducing sensitivity 
to these effects. A refinement to this is to restrict the off¬ 
set measurement to only a shorter, overlapped span of data 
between the two instruments, measure it in a separate fit pro¬ 
cedure, and hold the resulting value fixed in the main timing 
model fit. This will reduce covariance with long-term effects, 
but raises concern that the effect of the offset fit may not be 
fully accounted for in the other model parameter uncertain¬ 
ties. In both cases, the precision of the measurement is limited 
by the relevant TOA uncertainties. It is sometimes possible to 
transfer an offset measurement done using one bright pulsar to 
other sources, although there may then be concerns about po¬ 
tential systematics, for example due to pulse shape evolution 
with frequency, calibration inaccuracies, or different instru¬ 
ment configurations used for different sources. 

An alternate approach was used recently by Manchester 
et al. (2013), wherein a locally-generated pulsed signal was 
injected into the common signal path and used to recover the 
offsets between different systems with much higher precision 
than could be achieved using astronomical signals. This miti¬ 
gates all of the problems of TOA-based approaches described 
above. The only drawbacks are that it requires additional 
special-purpose hardware be built and installed at the obser¬ 
vatory, and that it can not be applied retroactively - once an 
instrument has been decommissioned it is no longer possible 
to perform this measurement. In contrast, simultaneous (or at 
least contemporaneous) pulsar data is often still available in 
archival data sets long after the relevant instruments are gone. 

We have developed a new method that addresses many of 
the shortcomings described above. This is based on the fact 
that for observations that are simultaneous in both time and 
frequency (i.e., where a single signal is split and fed into 
multiple backend systems), both instruments see not only the 
same pulsar signal, but also the same system noise. This cor¬ 
related noise can be used to recover a time offset with much 
higher precision than is possible from TOA-based methods. 
TOA determination can be viewed as a matched filtering pro¬ 
cess designed to recover the template profile shape. By con¬ 
struction, this filters out a large fraction of the noise that could 


Table 4 

Measured instrumental offsets. 0 


Receiver 

system 

Cross-corr 
offset (ns) 

J1713+0747 
JUMP (ns) 

J1909-3744 

JUMP (ns) 

Arecibo 327 

785(19) 

- 

- 

Arecibo 430 

789(5) 

- 

- 

Arecibo L-wide 

839(3) 

820(75) 

- 

Arecibo S-wide 

846(6) 

885(82) 

- 

GBT Rcvr_800 

897(8) 

951(124) 

936(42) 

GBT Rcvrl 2 

693(3) 

599(86) 

651(55) 


°Numbers in parentheses are uncertainties in the last digit quoted. 


J2145-0750 
J2010-1323 
B1937+21 
J1918-0642 
J1909-3744 
J1744-1134 
J1713+0747 
J1643-1224 
J1614-2230 
J1600-3053 
J1455-3330 
J1024-0719 
J1012+5307 
J0613-0200 


600 700 800 900 1000 1100 



Offset (ns) 


J2317+1439 
J2214+3000 
J2043+1711 
J2017+0603 
B1953+29 
J1944+0907 
B1937+21 
J1923+2515 
J1910+1256 
J1903+0327 
B1855+09 
J1853+1303 
J1741+1351 
J1738+0333 
J1713+0747 
J1640+2224 
J0030+0451 
J0023+0923 


i ru i 

MM S-wide 
MM 327 MHz 
HH 430 MHz - 
MM L-wide 


_... 1_ u _ _ 


1— 11 

. u+i .... 






.... J . J . . . .... 



r“1 , 

.LJ. J . . - I. 



.. | _ 




■ 




.. . . J . . . 




.. .hi _ 




in 

. L_f . . . .u. - - - 




. »■ ■ _ 




rrr—r 




._ 

. . . 1— ILJ _ 








.. . . 

rri 

.. .u .... u ... 




n ■ n 

. u ■.... 1 . 




n. i 



_i__]__i_i_ 


600 


700 


800 


900 


1000 


1100 


Offset (ns) 


Figure 4. Measured instrumental offsets versus pulsar for the GBT (top) and 
Arecibo (bottom). 


otherwise be used for an offset measurement. In our method, 
rather than cross-correlating measured pulse profiles with a 
noise-free template, pairs of simultaneous pulse profiles from 
each instrument are cross-correlated with each other. In con¬ 
trast with TOA measurement described in §3, it is advanta¬ 
geous to average the profiles as little as possible prior to this 
step, to preserve more (correlated) noise. With the profile 
data at its original time and frequency resolution, we compute 
cross-correlations between all pairs of profiles that overlap in 
both time and frequency. The cross-correlations from all si¬ 
multaneous profile pairs for a given pulsar, instrument setup, 
and epoch, are then averaged together (with a weight propor¬ 
tional to the amount of time-frequency overlap) to form a final 
correlation function. The lag for which this is maximized re¬ 
sults in an offset estimate for that portion of the data. We use 
this set of individual offset measurements to determine an un¬ 
certainty on the average offset, and look for systematic trends 
as a function of pulsar, time, or instrument setup. 
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For this work, we are interested in the offset between 
GASP and GUPPI at the GBT and the offset between ASP 
and PUPPI at Arecibo. We analyzed all available simulta¬ 
neously collected profiles with overlapping frequency bands 
from these pairs of instruments. The results of this offset 
analysis are shown in Table 4 and Figure 4. After account¬ 
ing for all a priori predictable latencies in the backend sys¬ 
tems, there remains ^700-900 ns additional offset between 
ASP/GASP and GUPPI/PUPPI that varies only as a function 
of signal path and instrument bandwidth. Table 4 lists the 
average value of all data available for each receiver system, 
while Figure 4 shows the same data averaged separately for 
each pulsar. At Arecibo, we obtain consistent results for the 
327 and 430 MHz setups, as well as for the L-wide and S- 
wide setups. This result is expected since these pairs of se¬ 
tups share common analog signal paths and PUPPI bandwidth 
(100 MHz and 800 MHz respectively). At Green Bank, the 
820 MHz and 1.4 GHz receivers have different signal paths 
to GUPPI, and the instrument is run at different bandwidth 
(200 MHz and 800 MHz). It is likely that both of these factors 
contribute to the observed offset difference. The sign of these 
values is such that ASP/GASP pulses arrive later - the off¬ 
sets must be subtracted from these TOAs to align them with 
the GUPPI/PUPPI data. As can be seen in Figure 4, consis¬ 
tent values are obtained from all pulsars for a given receiver 
system, therefore in our timing analysis we have applied the 
overall-average (Table 4) values to the TOAs. In our data set 
this is provided via a time offset (“-to”) flag on each TOA line. 

As a check on these results, we performed a standard timing 
analysis on the overlapping TOAs of two pulsars, in each case 
fitting for an offset between the TOAs from the two different 
instruments used as part of the timing solution. The results are 
shown in Table 4, where they are labeled JUMP (the tempo 
parameter used for this offset measurement). These values 
illustrate that our noise correlation provides both a consistent 
and much more precise result. For most other pulsars, TOA 
uncertainties are larger than for the pulsars used in Table 4, 
hence the JUMP uncertainties are larger as well. 


amplitude a and pulse phase shift 0 as 


X 2 (a,^) 


E 


\d k -at k e- 2 * ik *\ 2 


k 

D 2 +a 2 T 2 — 2aCdt (<j>) 


(Bl) 


where these terms come from the discrete Fourier transform 
of the measured pulse profile (< 4 ) and template profile ( 4 ): 


£ 2 =E l^l 2 

(B2) 

(N 

III 

<N 

Ex 

(B3) 

= t* k e 2 ^. 

(B4) 


k 


Here a 2 is the noise level in each bin of dk, and the sum is 
over pulse harmonics, not including the constant (DC) term. 
All pulse phase information is contained in the term, il¬ 
lustrating why TOA determination is sometimes described as 
a cross-correlation between the data and template profiles - 
the minimum x 2 is always achieved at the phase shift giving 
maximum cross-correlation. With the assumption of additive 
Gaussian noise (implicit in a x 2 fit), the TOA likelihood func¬ 
tion is 


p(d\a , 0 ) oc e 2 X = exp 


(2 aC dl (4>)-D 2 -a 2 T 2 \ 

V ^ y 


(B5) 


and with the use of uniform priors on a and 0 , the pos¬ 
terior distribution is simply proportional to the likelihood, 
p(a , 0| J) oc p(d\a, 0). For TOA determination, a is a nuisance 
parameter, which can be analytically marginalized over in the 
above expression to get the posterior 0 distribution 

p((/)\d) (X exp . (B 6 ) 


B. TOAS IN THE LOW-S/N LIMIT 

In the very low signal-to-noise ratio regime, the standard 
template matching procedure breaks down, producing under¬ 
estimated TOA uncertainties. In addition, the distribution 
of TOAs in this regime becomes significantly non-Gaussian. 
Here we derive expressions for the expected TOA probability 
distribution, and motivate our choice of S/N cutoff for TOAs. 
The use of a S/N or TOA uncertainty cutoff, or simply “by¬ 
eye” removal of outlier residuals, is often done in pulsar tim¬ 
ing analyses. The discussion in this section provides a some¬ 
what more rigorous and quantitative justification for this prac¬ 
tice. The behavior of TOA uncertainties in the low-S/N limit 
was previously explored empirically by Hotan et al. (2005) 
using simulated data, who reach similar conclusions to what 
we present here. 

We follow the standard Fourier-domain least squares TOA 
determination approach of Taylor (1992) (see also Demorest 
2007), writing the expression for x 2 as a function of fitted 


By making the substitution dk atk - i.e., the data profile 
is simply a scaled copy of the template - we can explore the 
expected shape of these distributions independent of any par¬ 
ticular (noisy) data realization. In this case Eqn. B 6 becomes 


/?( 0 ) oc exp 


/^Cg(0)\ 

V 2 T 4 ) ’ 


(B7) 


where S = aT/cr defines the signal-to-noise ratio of the data, 
and C tt is the template profile’s autocorrelation (with normal¬ 
ization C tt ( 0) = T 2 ). For non-detections (S —>> 0 limit), /?(0) 
becomes a uniform distribution across one turn of phase. For 
high-S/N detections (S > 10), p(0) becomes extremely well 
approximated by a Gaussian, with standard deviation given 
by the usual template-matching TOA uncertainty formula 

c r 4> = S~ l T (C'/i (0)) _1/2 . (B 8 ) 

In the low-S/N regime between these two limits, the standard 
uncertainty formula underestimates the true spread of TOA 
values and signifcant non-Gaussianity is present. We illus¬ 
trate this in Figures 5 and 6 using data from PSR J1455-3330. 
This pulsar provides a clear demonstration of this effect, be¬ 
cause its wide scintillation bandwidth and moderate average 
flux result in profiles with a large range of S/N values in our 
dataset. If included in a standard x 2 -based timing model fit, 
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Figure 5. Upper panel: 1.4 GHz template profile for J1455-3330. Lower 
panel: Expected pulse phase shift distributions for several values of S/N ratio, 
from Eqn. B7. For plot clarity these are normalized to 1 at </> = 0, rather than 
to constant integrated area. This shows the evolution of the distribution from 
nearly Gaussian at higher S/N ratio (S = 5) to clearly non-Gaussian at low-S/N 
OS = 2). 


the low-S/N points appear as outliers and have a dispropor¬ 
tionately large impact on the results. While it would be pos¬ 
sible to mitigate this by a modification of the TOA uncertain¬ 
ties or use of a timing model likelihood based on Eqn. B6, in 
our data set the amount of additional information gained from 
these data points is likely to be marginal at best. Instead, for 
the work presented in this paper we have simply removed all 
TOAs with S < 8. 



Timing residual / TOA uncertainty 


Figure 6. Signal-to-noise ratio S versus normalized timing residual (residual 
divided by TOA uncertainty) for J1455-3330 L-band data. All data points 
below the solid line at S = 8 were removed from the timing analysis. 


matrix that has columns of alternating sine and cosine func¬ 
tions for frequencies in the range [1 /T, n mo fe/T ]; T is the total 
observation time span, A/ = 1 /T, and n mo d e is the number of 
frequencies included in the sum. The vector a gives the am¬ 
plitudes of the Fourier basis functions (see Lentati et al. 2014; 
Arzoumanian et al. 2014, for more details). The term f/j de¬ 
scribes noise that is uncorrelated in time but completely cor¬ 
related between TOAs obtained at different frequencies mea¬ 
sured simultaneously. This term could be due to pulse phase 
jitter but could also have other components not due to jitter. 
This term characterizes noise that is completely correlated for 
all TOAs in a given time bin but completely uncorrelated be¬ 
tween time bins. The matrix U is an A/joa x N t b matrix that 
maps TOAs to a given time bin and j is the amplitude of the 
short time-scale fluctuations. Finally the last term n describes 
a Gaussian white noise process that characterizes time-, and 
frequency-independent random noise left in the data. 

Since the white noise is modeled as Gaussian, the likeli¬ 
hood function for the noise is given by 


where 


exp (-^n t N ] n 
Vdet(27rA5~’ 


Nij^EpWij + QlStj), 


(C2) 

(C3) 


C. DERIVATION OF NOISE MODEE LIKELIHOOD 

We begin by forming a set of residuals via the standard 
weighted least squares fitting routine. An A/toa length vec¬ 
tor of residuals can be modeled mathematically as the sum of 
several deterministic and stochastic sources as follows 

£t = M€+Ea + £/j + n. (Cl) 

The first term on the right hand side (Me) describes small de¬ 
terministic trends due to timing model subtraction. Here M is 
the timing model design matrix and e is a vector of small tim¬ 
ing model parameter offsets. Next, the term F a models the red 
noise via a Fourier decomposition 10 - F is the Fourier design 


is an A/joa x A TO a matrix with E k and Q k corresponding to 
tempo and tempo2’s EFAC and EQUAD parameters for 
each observing backend, respectively, W = diagjof} is a di¬ 
agonal matrix of TOA uncertainties, and dij is the Kronecker 
delta function. The notation is such that the matrix elements 
(/, j) apply to those TOAs corresponding to the backend ob¬ 
serving system labeled by k. We can now write the likelihood 
function of the residuals as 


p(<5t|e,a,j,0) = 


exp (-\r T N *r) 


\J det(27rV) ’ 

where <fi denotes the E k and Q k parameters and 


(C4) 


r = St-Me-FsL-Uy (C5) 


10 The Fourier basis was chosen to improve computational efficiency; it is We now wish to impose prior distributions on our short 

not a requirement of this noise modeling method. timescale correlated noise and red noise. Both can be mod- 
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eled as Gaussian processes by imposing the following priors The marginalized posterior distribution is then 


„/:m_ ‘j) 

v i 

(C6) 

P( a| A) = exp <±!E! a) , 

v/det(27 Tip) 

(C7) 


where = J%5ij is an N tb x N tb matrix with diagonal ele¬ 
ments, and Jj describes the variance of the jitter-like corre¬ 
lated noise for each observing backend; it is also referred to 
as the ECORR parameter in tempo and tempo2. Further¬ 
more (fij = diag{10 p "} is an 2ft mode x 2ft mode matrix describ¬ 
ing the variance of the red noise Fourier coefficients at each 
frequency. In this framework, the coefficients of the (^-matrix 
are related to the power spectral density evaluated at a given 
frequency. In principle we could use the power spectrum co¬ 
efficients, 10 p % themselves as free parameters but in this work 
we parameterize them via a power law 

<p«*w A. = J_ A 2 YAV red (C8) 

Apan VJyr/ 

where r span is the total observation time, A red is the amplitude 
of the red noise in /is yr 1 / 2 , 7 red is the spectral index of the 
red noise, / yr is the reference frequency of 1 yr -1 , and f n is 
the ftth Fourier frequency assuming Nyquist sampling. We 
see that the prior distributions on jitter-like correlated noise 
and red noise are themselves parameterized by some combi¬ 
nation of hyper-parameters. We can write down the posterior 
distribution for the residuals 

P(e,a,j,(f)\5t) oc p(5t\e,a,j,4>)p(j\J k )p(a\p n ). (C9) 

For the purposes of estimating the underlying noise charac¬ 
teristics of our data set, the parameters e, j, and a are nuisance 
parameters that we wish to marginalize over. This can be done 
in a sequential fashion as was presented in Arzoumanian et al. 
(2014), but here we take a different approach. Notice that all 
timing parameters are linear and can be described with Gaus¬ 
sian prior distributions 11 . We can thus define a combined op¬ 
erator matrix and amplitude vector 


p{(j) |£t) oc / 

J-GO 

oo 

dbp(.5t\b)p(b\ct>)p(<fi) (C14) 

OO 

exp [-iO^AT^t-d^ir'd)] 

_ V / ( 27 r) A, TOA-dimb det(Af) d e t(fi)det(E) ’ 

where 

d = T T N~ l 5t (Cl 5) 

E = (B~ l + T t N~ x T). (Cl 6) 

The maximum likelihood values of b and their uncertainties 
can be found as 

b = E _1 d (C17) 

cov(b) = E _1 . (Cl 8) 

This scheme has the advantage of being computationally ef¬ 
ficient in that it bypasses 0(Nj OA ) matrix operations via rank 
reduced matrices (van Haasteren & Vallisneri 2015) resulting 
in a likelihood evaluation that instead scales as 0(A par ), where 
A/par is the sum of the number of timing parameters, red noise 
sample frequencies, and observing time bins. For the largest 
datasets the computational speed up is a factor of ^ 10 3 . 

For a given set of hyper-parameters, this allows us to deter¬ 
mine the maximum likelihood timing model parameters and 
the maximum likelihood red noise realization present in the 
data via the equivalent of a generalized least squares fit. We 
can also evaluate the posterior of the hyper-parameters 0 and 
thus find the maximum likelihood noise parameters including 
the EFAC, EQUAD, ECORR, red noise amplitude A red and 
spectral index 7 red . The posterior distributions of the noise 
parameters are sampled using a Markov Chain Monte Carlo 
process in which we sample some parameters in log 10 space 
and limit them to log 10 *4 £ [-8.5,-4], with A in units of sec¬ 
onds, log 10 A red G [-7.5,1.5] where A red is in units of /is yr 1 / 2 , 
and 7 red G [0,7]. 


T = [M F U] , 



(CIO) 


with prior distribution 


P( b|0) 


exp (-|b t B 1 b) 
V"det(27r B) 


(Cll) 


and covariance matrix defined in terms of the block matrix 


B = 


oo 0 0 

0 ip 0 

0 0 J 


(Cl 2) 


where oo is a diagonal matrix of infinities to describe a uni¬ 
form prior on e. The resulting likelihood function is then 


p(<5t|b) = 


exp [-±(5t-Tb) T N-\5t~Tb)] 
Vdet(27rA0 


(Cl 3) 


11 We use uniform priors on the timing model parameter offsets, e but this 
is the same as a Gaussian prior with infinite variance. Technically this prior is 
not normalizable, but since we are interested in parameter estimation and not 
Bayesian model selection here, this non-normalizable prior is not a problem. 


D. DAIFY AVERAGED RESIDUAFS 


For modem wide-band timing campaigns using multi¬ 
channel TOAs it becomes useful to visually inspect timing 
residuals that have been averaged in order to look for long 
term trends or biases. Here we derive a robust weighted av¬ 
erage that will fully account for short timescale correlations 
introduced by the ECORR in our noise models. This is im¬ 
portant since ECORR is meant to model pulse phase jitter, 
thus when constructing daily averaged residuals, one must in¬ 
clude this effect as it results in larger averaged uncertainties 
on the averaged residuals. In essence this allows for a way to 
visually determine which pulsars may be dominated by pulse 
phase jitter. 

We begin the derivation by introducing the probability dis¬ 
tribution of the group of residuals that belong to time bin 12 
k, 


p(St k \St k ) = 


exp[4(r)t, - OShfCi 1 ((51, - OSt k )] 
det(Q) 


(Dl) 


12 In this work, we have used time bins of size 1 second, thus are only 
averaging sets of multi-channel residuals measured simultaneously. 
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where Stk, Stk, and Q are the residuals in time bin k , the mean 
residual in time bin k , and the covariance matrix of the residu¬ 
als in time bin k , respectively. Here, O is the design matrix for 
the mean which in this case is a vector of ones of length Nk, 
where Nk is the number of residuals in simultaneous time bin 
k. In an identical manner as Appendix C we can determine the 
maximum likelihood estimator and uncertainty for the mean 
of the probability distribution function (i.e. the daily averaged 
residual) 


Stf h = (0 r C- k '0)-'0 r C- k 'St k (D2) 

cr 2 k = (0 T C- k '0)-\ (D3) 

where is the weighted uncertainty on the daily averaged 
residual. Note that if Q is diagonal with elements correspond¬ 
ing to the TOA uncertainties then we obtain our usual expres¬ 
sion for the weighted mean and standard deviation 








i,k 



al=i^nT 2 


(D4) 

(D5) 


where 07 ^ is the TOA uncertainty for the i TOA in simulta¬ 
neous time bin k. We note that the ECORR will add to the 
off-diagonal components of Q and can have a large impact 
depending on the relative strength of ECORR compared to 
the radiometer noise component. 
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Figure 7. Timing summary for PSR J0023+0923. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 8. Timing summary for PSR J0030+0451. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 9. Timing summary for PSR J0340+4130. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 10. Timing summary for PSR J0613-0200. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 11. Timing summary for PSR J0645+5158. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 


J0931-1902 


GASP GUPPI 



Figure 12. Timing summary for PSR J0931-1902. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 13. Timing summary for PSR J1012+5307. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 14. Timing summary for PSR J1024-0719. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 















































































NANOGrav Nine-year Data Set 


23 


J1455-3330 


GASP 


GUPPI 


I 80 


£ -80 


<D 1 - 1 

a 

u ^ 

9 *6 

< s 

& 


KS S 

X o 

S a 

Q 


1.5 

o.o k 


t - 1 - 5 

^ -3.0 


“I-1-r- 


T 


-- 


“i-i-r 


T 


i i- r~ 


| - | - 





H- h 


H-1-h 


o 7* — — ^\\i\ -I— 4 ;t |-j |t -j*H ilj - Hr♦'* jf* — : 

-20 r 1 1 I 1 t 

!-1—M—H-1-1-1 H-1-1 H--M—I-I— 1 —I-1-1-1-1-1-! 


il I l 


H-h 


H-|H 



2006 


2008 


2010 
Date [yr] 


2012 


2014 


Figure 15. Timing summary for PSR J1455-3330. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 16. Timing summary for PSR J1600-3053. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 17. Timing summary for PSR J1614-2230. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 18. Timing summary for PSR J1640+2224. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 19. Timing summary for PSR J1643-1224. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 20. Timing summary for PSR J1713+0747. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 21. Timing summary for PSR J1738+0333. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 22. Timing summary for PSR J1741+1351. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 




























NANOGrav Nine-year Data Set 


27 


J1744-1134 


GASP 


GUPPI 


a 
P 
T3 
c ~ri 
<D 
Pi 


<D 1 — 1 

o3 P 
5h ^ 

£ ^ 
< s 

Pi 


I 

^ s 

X o 

2 & 

Qcn 

o 


15 

0 

-15 

3 

0 

-3 

0.8 

0.4 

0.0 

-0.4 

- 0.8 


:-'-'- 1 - 

- • 

_ 

i i i • | 1 1 | 

^—-—-—r— 1 —»—r—i ; 

liiaiilihiiJiiiUnfc iiMiii 

t# * *• * — 

J 1 X I X X L_ 1 : 

1 

r 

i 

i 

_—u— _ 

% 

Id-Hty- - ■< 

i i i 

^ jt-j -*|tj 

j_ 1 _ 1 

~±~ 

- 

1 

J_I_L. 

**V***k***rt -: 

i_i_i_i_I 

1 1 

f 1 * 

:_i_i_ 1 _ 

1 1 1 

I 

1 1 1 

I 

J_i_i_i_L_ 

1 i i i i i i : 

1 

1 _ “ 

. i -jt 

| —MT-n- 


2006 


2008 


2010 
Date [yr] 


2012 


2014 


Figure 23. Timing summary for PSR J1744-1134. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 24. Timing summary for PSR J1747-4036. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 25. Timing summary for PSR J1832-0836. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 26. Timing summary for PSR J1853+1303. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 27. Timing summary for PSR B1855+09. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 28. Timing summary for PSR J1903+0327. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 29. Timing summary for PSR J1909-3744. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 









































NANOGrav Nine-year Data Set 


31 


J1910+1256 


ASP 


PUPPI 


cd 

P 

T3 

<U 

Pi 


<D 1 — 1 

2 § 

< SS 


<D '—' 

P ^ 

O cd 
.tn =3 

^ <D 

P< 


I 

K> S 

X o 

S a 

Qcn 

o 


-15 

6 

0 

-6 

3 

0 

-3 

2 

0 

-2 

-4 


t -****-•♦*£-- * Ll ^ 

— , 

• 

:..1...1..,1...1 

1— T - 1 - r ~^“l -: 

H#niSH 

1 •••?•* 

1 1 I* 1 1 1 

: i i 

-k- 

- V 

:i_i 

1 1 1 1 

1 1 1 1 1 1 1 

1 1 11 

i 1 i i i 1 

| 1 I I | : 

1 -= 

■t i 

1 i i i 1 

1 


WlM-HPf-l- 

fl 1 1 1 1 1 

1 1 1 'll: 

1 ' T 

j|t|-))44|##J(K- - 
1 - 

l i i i i i i i 

i 3 

L I T 

__1_I_1 1 '_' 1_1_1___ 

itfaisi M l 

j_1_i_i_i_1_ 

- & - 

■ 1 ' i i 1 i i i 1 i i i 1 i i i 


2006 


2008 


2010 
Date [yr] 


2012 


2014 


Figure 30. Timing summary for PSR J1910+1256. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 31. Timing summary for PSR J1918-0642. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 32. Timing summary for PSR J1923+2515. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 33. Timing summary for PSR B1937+21. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 34. Timing summary for PSR J1944+0907. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 35. Timing summary for PSR J1949+3106. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 36. Timing summary for PSR B1953+29. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 37. Timing summary for PSR J2010-1323. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 38. Timing summary for PSR J2017+0603. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 39. Timing summary for PSR J2043+1711. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 40. Timing summary for PSR J2145-0750. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the top 
panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 41. Timing summary for PSR J2214+3000. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 42. Timing summary for PSR J2302+4442. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 
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Figure 43. Timing summary for PSR J2317+1439. Colors are: Blue: 1.4 GHz, Purple: 2.3 GHz, Green: 820 MHz, Orange: 430 MHz, Red: 327 MHz. In the 
top panel, individual points are semi-transparent; darker regions arise from the overlap of many points. 





























